A self-organized critical model and multifractal analysis for earthquakes in Central Alborz, Iran

This paper is devoted to a phenomenological study of the earthquakes in central Alborz, Iran. Using three observational quantities, namely the weight function, the quality factor, and the velocity model in this region, we develop a modified dissipative sandpile model which captures the main features of the system, especially the average activity field over the region of study. The model is based on external stimuli, the location of which is chosen (I) randomly, (II) on the faults, (III) on the low active points, (IV) on the moderately active points, and (V) on the highly active points in the region. We uncover some universal behaviors depending slightly on the method of external stimuli. A multi-fractal detrended fluctuation analysis is exploited to extract the spectrum of the Hurst exponent of the time series obtained by each of these schemes. Although the average Hurst exponent depends slightly on the method of stimuli, we numerically show that in all cases it is lower than 0.5, reflecting the anti-correlated nature of the system. The lowest average Hurst exponent is found to be associated with the case (V), in such a way that the more active the stimulated sites are, the lower the average Hurst exponent is obtained, i.e. the large earthquakes are more anticorrelated. Moreover, we find that the activity field achieved in this study provide information about the depth and topography of the basement, and also the area that can potentially be the location of the future large events. We successfully determine a high activity zone on the Mosha Fault, where the mainshock occurred on May 7th, 2020 (M\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_W$$\end{document}W 4.9).

It is widely believed that the earthquake is a self-organized critical system. When an earthquake occurs and a fault slips, it causes an excess in the tension of the neighboring regions. The motion of the neighboring regions depends on their local strain, so that if it exceeds a threshold, then it slips. This dynamics is reminiscent of the sandpile dynamics, first invented by Bak-Tang-Wiesenfeld (BTW) model, where the energy spreads the system based on similar dynamic rules: when the local energy exceeds a threshold, the site topples, rising the energy of the neighboring sites by one unit. BTW and some other variants 1,2 are too unrealistic to give an acceptable description of the earthquake, e.g. it cannot explain 1 f noise that is seen in reality 3 . For more realistic situations, one needs more detailed model paying attention to the structure of the earth and also the information on how the seismic activities affect each other, which translates to how a perturbation propagates from region to region. The latter is crucial and pretty complicated since it depends on the material content of the Earth's interior where the signal propagates. Recently, using a virtual seismometers, the correlations between seismic activities were incorporated in the model, based on which a complex network was designed on top of which sandpile dynamics were implemented 4 . The virtual seismometer can provide inter-event empirical Green's function in the Earth's interior. Therefore, using unconventional form of seismic interferometry, one earthquake beneath the Earth's surface can turn to a receiver where recorded another event waveform 5 . One may think of this problem from another point of view: the cross-correlation between previously happened earthquakes gives us a set of valuable information about the structure of the earth and the signal propagation in the region. For example, suppose that we have a map (the place and the magnitude, M) of previous earthquakes in the region of interest, along with the corresponding time series. Cross-correlation between the events gives us a criteria of how events are related and to which extend they are correlated.

Observational data, modified sandpile model
A network with evenly grid space was used in this study, which has been achieved by processing the waveforms of the occurred earthquakes in the Central Alborz region. Calculating the weight function between cells involves the interferometry of the recorded waveforms, attenuation, and seismic velocity model of the study area. The workflow of the paper is shown in Fig. 1. As is evident in this figure, we follow three objectives in this work: • 1: Using the velocity model, Q-factor, and the weighted matrix between cells (calculated by signals-crosscorrelation procedure), and the phenomenological parameters, we construct an effective correlated lattice. • 2: We developed an effective dynamical (sandpile) model for the earthquake. Earthquake dataset. We processed all tectonically earthquakes that occurred in the Central Alborz with M ≥ 2.5 (upon which a higher concern in the community is), between 2006 and 2021 except an earthquake occurred on May 7th, 2020 (M W 4.9). The reason for considering only this range is that the earthquakes with smaller magnitudes are often (but not necessarily) aftershocks, the location of which is almost a priori known. The other reason is that we do not have access to the catalog of these small events. These events were recorded by 48 seismic stations operated with three permanent seismic networks that included the following: (I) Ira-  Figure 1. The work flow of the paper. We start from the data set obtained by the velocity model, the Q-factor, and the weighted matrix between cells (calculated by signals-cross-correlation procedure). These quantities are employed as the input parameters of a modified sandpile model. At the end, the time series of the sandpile model is studied using the multifractal analysis of detrended time series. www.nature.com/scientificreports/ malizations to suppress the influence of instrument irregularities, human activities, and source time functions with nonuniform energy. The maximum amplitude of the envelope function within the expected signal window (1.5-3.5 kms −1 ) to the root-mean-square, rms, is the SNR definition in this study 32 . The amplitudes of the waveform parts beyond the expected signal window were transferred to zero. Next, we took the waveform part from the origin time, t o , to the end of the Rayleigh coda wave (1.5 kms −1 ) for the signals-cross-correlation procedure. Then, we applied the signal-cross-correlation (hereafter SCC) operator on the prepared waveforms of a pairevent which are recorded by a common station in the alignment of the inter-event line. The total number of inter-event raypaths depicts in Fig. 2b which is 47,864 paths. The SCC can be summarized where T m , ∂ , Ŵ , r, u, and T are the moment tensor solution, spatial gradient, inter-event Green's function, coordinate vector, displacement and traction, respectively. Moreover, the indices 1 , 2, and S represent event1, event2, and Earth's surface, respectively. Since we have taken into account only vertical components of signals, the full moment tensor solutions reduce from nine components , which is a constant value. The effect of these constants ( T m ZZ 1 and T m ZZ 2 ) have been removed using a simple time-domain normalization operator in the data preparation step. Then, the maximum amplitude of the absolute SCC function, SCC A max , was extracted and attributed to the corresponding cells along the inter-event raypath. Given the limitations and non-uniform distribution in raypath coverage, we used a grid base tomography procedure to obtain SCC weight for all evenly grid space for the area under study. For tomography, the observed data, d obs , is t SCC A max (the arrival time of SCC A max ) and model, m cal , could be the SCC A max . Therefore, the relation between data and model is The calculated data, d cal , can be calculated using Fast Marching Method (FMM; 33,34 ) grid base algorithm. When G, Green's function, is a linear or near-linear function this formula can be An iterative linearized damped-least squares inversion procedure can apply to minimize the observed and calculated data misfit 15,35 . Using Gauss-Newton gradient method, the relation (3) can be where ǫ and η are the damping, and smoothing regularization parameters, and also C −1 d and C −1 m remark the data, and model errors, respectively, and Tr represent the transpose operator. To solve equation 4, the area under study was meshed with the grid size 14 km × 14 km, and the average of A max applied as initial input model, m 0 . Also, the regularization parameters, ǫ , and η , were obtained by standard the L-curve by considering a trade-off between data misfit and model roughness. To stabilize the result, we used those observed data, A max , with residuals less than two standard deviations (2σ A max ) in the inversion procedure. Figure 3a shows the obtained SCC weight map in the study area. The damping value, ǫ , was fixed to 255, and the smoothing parameter, η , was 950.
The first arrival P-wave traveltimes can be used to obtain the 3D crustal velocity structure of the Central Alborz region. This 3D model can provide an insight into the detailed crustal velocity structure of the Central Alborz to better understand the fine-scale tectonics and seismic data transfer speed. Based on the event's epicenter distribution and our study area, the 3D velocity model obtained by Afra et al. 36 can be an appropriate model. This model has been calculated using an iterative linearized, damped least-squares widely used inversion code SIMULPS 37 . Because the recovered anomalies of the model have been confirmed by the two different (including seismology and gravity) geophysical methods. Moreover, the reliability of this velocity model was performed by different resolution tests which are including the checkerboard resolution (both dense and sparse for checking lateral and smearing resolutions, respectively), input initial model uncertainties, events location uncertainties, and Resolution Diagonal Element 36 . In this study, the region with the fairly resolution is surrounded by a black thick border according to 36 . Fig. 3c represents this 3D velocity model from subsurface to depth of 30 km. It's worthy to note that the 2D SCC weight was also parameterized as the depth of the velocity model as shown in Fig. 3b.
The quality factor (Q-factor) generally refers to the attenuation of the amplitude due to elastic and inelastic properties of the medium. The coda and body wave Q-factors are a powerful tool to study thermal, compositional, and deformational characteristics of Earth's interior 38 in seismology. The body or surface waves are attenuated with rates greater than the calculated rates for geometrical spreading. The inverse of the Q-factor is composed of two factors, namely the scattering Q −1 Sc , and the attenuation along the propagation path of a seismic wave (due to the inhomogeneities within the earth and the geometrical spreading) 29 obtained the Q-factor model for the Central Alborz. This model has been calculated by the Lg coda method using 1020 waveforms of the vertical component of 205 earthquakes with 3.5 ≤ M L ≤ 6.5 recorded by 35 short-period seismic stations between 2000-2009. In this study, we used the Q-factor model calculated by Naghavi et al. 29 for predicting source-receiver attenuation as depicted in Fig. 3d. www.nature.com/scientificreports/

The dynamical model
Our dynamical model is a variant of dissipative continuous sandpile model 40,41 , that is implemented on top of a three-dimensional cubic lattice with coordination number z = 6 . The spatial extent of the region under study is the map presented in Figs. 2 and 3. We parameterized (meshed) the system so that the lattice points are fitted to existing data from which we picked the earth factors like the q-factor and the velocity filed. The resulting lattice consists of N 1 × N 2 = 100 × 70 nodes in each plate parallel to the earth surface, and totally N 3 = 13 horizontally plates are considered in the perpendicular direction as shown if Fig. 3, so that the lattice includes N = L 1 L 2 L 3 sites. To each site of the lattice i, we attribute three intrinsic quantities obtained from the real data of the earth explained in the previous section 29, 36 : the weight function ( W i ), the quality factor ( Q i ), and the velocity Model ( V i ). The weight and other functions have designed so that their values for the connection between two neighboring sites is The resulting weight field, the Q-factor, and the velocity field model are shown in Fig. 3b, c and d. Using the weight field we construct the lattice as represented in Fig. 3a, where the position of faults and the points at which earthquakes have taken place are shown. The threshold field to be used in the dynamics of the system is proportional to the field given in Fig. 3b, i.e. ǫ th i = W ij . Once the network is constructed, we define the following dynamics, which is based on the local energy/stress in each site i, denoted by ǫ i , taking values in the range [0, ǫ th i ≡ z j=1 W ij ] (the summation j is over neighbors of i), so that a local status of the system is identified by the set {ǫ i } N i=1 . The initial state of the system is chosen to be random with uniform distribution.
The dynamics of the system are defined by local relaxations generated by local slipping of the fault, i.e. distributing the stress excess through the neighboring regions. The rate of stress transfer is related to the Q-factor, velocity model, and the weight of the connections. The local stimulation of a region is external and is implemented on a randomly chosen site i via ǫ i → ǫ i + r , where r is a flat random number between 0 and 1), which favors a local slip of this site. This site is however static if its accumulated stress is lower than a threshold ǫ th i as a consequence of local static friction 42 . If ǫ i exceeds ǫ th i , then site i is called unstable, leading to a toppling process (local relaxation), during which ǫ i → ǫ i − � ij , where matrix is defined as where A 0 is an amplitude, f is frequency, r ij and t ij ≡ r ij V ij are the distance and the travel time between sites i and j. In our coarse grained model, we ignore the dependence on the frequency by setting f = const. . Moreover, to recover conservative dynamics in the limit Q → ∞ , we set A 0 = √ r ij , so that the inelastic attenuation factor reduces to A ij = exp −r ij /Q ij V ij . The toppling rule is schematically shown in Fig. 4. Also we have r ij = 4.52 km in the x 1 direction, r ij = 5.55 km in the x 2 direction, and r ij varies from 2 km to 4 km in the x 3 direction.
Burst dynamics is a popular property of sandpiles, manifested by prominent avalanches which occur with a low frequency, called sometimes rare events [43][44][45][46][47][48] . Avalanches in our model are defined as a chain of local relaxations (topplings) occurring as a consequence of an external stimulus. In fact, when one site becomes unstable by an external stimulus, it topples, and as a consequence, the neighboring site might become unstable and topple in their turn, so that a chain of relaxations take place up to a time where no further unstable site is found. The duration D and the size S of avalanche are defined as the lifetime and the total number of topplings in the avalanche respectively. Then another random site for stimulation is chosen and so on. The average local stress grows almost linearly with time until reaching a stationary state after which the amount of stress that leaves the system through the boundaries is statistically equal to the number of input stress, for a good review see 49 .

Measures and results
In this section, we present the results of the simulation of our modified sandpile model. The size of the lattice is fixed as explained in the previous section. We have tested five kinds of external stimuli, which are listed bellow: • I: completely random stimuli (the sites for external stimuli is completely random chosen),   www.nature.com/scientificreports/ The case I represents the case were the local stress excess (resulting to an earthquake) takes place in a completely random region, while the case II realizes the situations where the earthquake starts on the fault. The cases III, IV, and V capture the cases where its starts from more active sites which is classified to M ≥ 2.5 , M ≥ 3.0 and M ≥ 4.0 . For each case, we have generated over 10 6 samples, i.e. the avalanches in the stationary states. The activity field for 10 6 samples for the cases I, II and III are shown in Figs. 5, 6, and 7 respectively. By looking at Fig. 3a, we see a good correlation between the average activity field and the weight field. A much similarity is observed  www.nature.com/scientificreports/ between the case I and the weight field but as the pattern of stimuli changes (cases II, III, IV, and V), the activity field show different patterns. Especially for the fault stimuli (case II) the activity in the vicinity of the fault positions is much higher, which is rather expected. A more realistic situation is the case where the stimulation takes place in the vicinity of the active points (cases III, IV, and V), e.g. Fig. 7 since these points are more active with respect to the rest regions. Apart from the observational data (the points where the earthquake has taken place), this activity field predicts the activity of the other regions, and explicitly shows important points which can be potentially the starting point of the upcoming earthquakes. This shows that our model supports the location of the next earthquake, while the input parameters (cross-correlation weight, velocity model, Q-factor, etc.) of our processing did not include the information of this earthquake. Although we applied the random event location (results depicted in Fig. 5), the potential to produce an earthquake mostly can be appeared around Kuh-Sorkh and Parchin faults. By exciting cells on the known faults (red lines in Figs. 1 and 2), this potential mostly appears around Atari-Firuzkuh, western of Mosha, Kuh-Sorkh, and Parchin faults. This potential can become apparent for cells excited by events' location around Atari-Firuzkuh, Firuzkuh, and western of Mosha. A simple comparison for these activities (Fig. 5, 6, and 7) indicate the potential to produce an earthquake for superficial layers can be expected for micro-earthquakes (M ≤ 3.0) which is in agreement with seismicity in the study area. But, this potential can produce a larger earthquake at the greater depths (15 to 30 km) as experienced by historical earthquakes (see Fig. 2). Our results are consistent with previous studies (e.g. 36 ) that P-wave tomography and gravity inversion models reveal a region at depth ranges of 12.5-17.5 km (follow the orange contour in Fig. 8 between distances of 60-90 km) for the potential of producing an earthquake with M 6.5.
The left-lateral Mosha fault, which is the most important internal fault to the central Alborz Mountains, can release a part of energy in the amount of the mm yr −1 range-parallel strike-slip motion. Vertical intersection distributions of the final activity models for Mosha fault are shown in Fig. 8. This profile along ( 51.50°E, 35.87°N ) to ( 53.30°E, 35.60°N) is ∼ 165 km long. The random activity along this profile can clearly represent the upper crust (a thick layer up to 20 km depth) in which most earthquakes can occur. This thickness is in agreement with the report by Abbassi et al. 18 obtained by the employment of a local seismic network. For the fault activity, the interface topography between the upper and middle crust improves realistically, which is in agreement with the bottom of event depths as reported in 50 . Inspection of this intersection reveals that high activities seem to be around the Firuzkuh, Atari faults, and also Mosha-North Tehran fault intersection (M.N.T.I) in an expectation area (see 36   www.nature.com/scientificreports/ The results that are shown in the above figures are graphical demonstration of the situation that the region have. It is now worthy to quantify the universal behaviors of the model for the five cases that we introduced. In Fig. 9, the distribution functions (P) of the avalanche duration D and size S are exhibited. In Fig. 9a and b, we show the log-log plot of P(S) and P(D) respectively where the linear decrease in the signature of power-law decay P(x) ∝ x −τ x , x = S, D . Interestingly the exponents depend of the taken situation, see TABLE 1. Noting that for the three-dimensional Bak-Tang-Weisenfeld (BTW) sandpile model, we have τ S ≃ 4 3 44,47 , we see that for all cases, the exponents τ S and τ D coincide well with 3D sandpile. The scaling dimension γ SD defined by S ∝ D γ is almost robust, i.e. γ SD = 1.75 ± 0.008 , 1.73 ± 0.008 , 1.76 ± 0.009 , 1.74 ± 0.008 and 1.74 ± 0.008 for the cases I to V respectively.
To test the criticality of the system, one can use the branching ratio function defined by the conditional expected value b(x) ≡ E S t+1 x |S t = x , where E is expected value. For the critical systems lim x→0 b(x) = 1 or is in the vicinity of unity 4,51-54 . This function has been shown in Fig. 9d, for which 1 < lim x→0 b(x) < 1.03 , showing that although not exactly (super-critical), but we are in pretty in the vicinity of the critical region, which is confirmed by other power-law behaviors.

Multifractal detrended fluctuation analysis
Here we analyze the time series of our sandpile model and inspect the spectrum of the Hurst exponent. The method that we describe in this section is multifractal detrended fluctuation analysis (MDFA), a technique most commonly used in the analysis of self-affine fractal time series. For an uncorrelated time series the Hurst exponent is a single value H = 0.5 55,56 . For a mono-fractal time series, the Hurst exponent, H, is defined in terms of the asymptotic behaviour of the rescaled range. Consider a general time series {x(t)} n t=1 , and the cumulative deviate (profile) series The range is then defined by The standard deviation is also defined as For multi-fractal time series, different Hurst exponents govern the system in different size scales. The spectrum (importantly the average and the width of the spectrum) of the Hurst exponent plays the dominant role in identifying the properties of the time series of earthquake. The multifractal analysis (MFA) is an important tool in such time series (which are not described by a single Hurst exponent and single scaling relations). Apparently, the inspecting such systems using single scaling relations, gives, at best, the average Hurst exponent, and sometimes a wrong result. Even if clean scaling relations are observed using this method, it does not say anything about the other exponents available in the system for different size scales. MFA instead gives us the spectrum of the Hurst exponent for these time series. Therefore, in this section we asses the multifractal analysis, which studies the earthquake time series (which is multifractal) with a spectrum of the Hurst exponent (for the monofractal time series, the spectrum is peaked with a zero width). This spectrum lets us know the mean as well as the fluctuations of the Hurst exponent, helping us to distinguish the type of correlations (and anti-correlations) existing in the system.
To describe MFA, we define a variance for each of the segments v = 1, 2, 3....N s by 57 www.nature.com/scientificreports/ where Y (v) represents the mean of Y over the segment v. The q th moment is then obtained using The typical behavior of F q (s) is as follows 58 where h(q) is the corresponding exponent that is related to the Hurst exponent. To extract the spectrum of the Hurst exponent we use the standard multifractal analysis. Therefore, we use the generalized variance Eq. (10), but now with a new exponent where τ (q) represents the classical multifractal scaling exponent. This exponent is related to h(q) for stationary and normalized time series, The Legendre transform of the generalized scaling exponent τ (q) gives the multifractal function as follows where α = ∂τ (q) ∂q . Finally by a simple replacement one obtains which gives the spectrum of the Hurst exponent. The function F(S) is shown in Fig. 10 for the cases I till V and for various amounts of q , where a power-law behavior is evident in a large interval (nearly two decades). The exponents of these graphs are h(q) . If h(q) is the same for all q values, then we have mono-fractal with Hurst exponent h(q = 2) . For this case f (α) would be a peaked function around α = h(q = 2) with zero width. Figure 11 shows that this is not the case, and we are facing with a strong multifractal time series for all cases. Figure 11a shows the result for which we let the avalanches go beyond the almost square region identified in Fig. 2b (high-resolution region), i.e. the region with high resolution that we are more confident about the weight field that we obtained. Figure 11b shows the results for the case where we restrict the avalanches to the square. For both cases, we see that the width of f is pretty high, and the peak position varies with the method of stimulation. Even for completely random stimulation, the peak is around ᾱ = 0.37 (0.36) for the avalanches in the whole space (inside the high-resolution box), both being lower than 0.5, showing that the system is anticorrelated. The exponents are shown in the Table 1. Using the fact that the exponent of autocorrelation function ξ is related to the Hurst exponent like ξ = 2 − 2H 59 this result uncovers that the time series are anti-correlated, meaning that a large event is often followed by a small event and vice versa. This effect magnifies when the stimulation is more selective, i.e. for the fault stimuli, it is ᾱ = 0.36 , and for highly active stimuli case ( V ) it is even smaller, ᾱ = 0.33 for the case where avalanches are allowed to go through all the space. The main results of our analysis are gathered in Fig. 12a and b (in whole space and limited to the high-resolution box respectively), where the relation between various exponents are shown in terms of the strategy of stimuli. The decrease of the Hurst exponent (green diamond symbols) with the stimuli strategy is evident in these graphs, reflecting the fact that the time series becomes more and more anti-correlated. Along with this, the exponents τ S and τ D decrease, but γ SD are almost robust.

Conclusion
In this study, we focused on the earthquakes in central Alborz, Iran. In the first part of the paper, we explored the properties of the earth in the region under study, as well as the rate of earthquakes. It helped us to construct a modified sandpile model which is much similar to the continuous dissipative sandpile model in which the energy dissipation is related to the quality factor and the velocity model of the earth. The weight function which was obtained using the signals-cross-correlation of the real seismic activities was used to estimate the weight field which was employed for distributing the energy to the neighboring sites in each toppling. Our model is based on external stimuli, the location of which can be (I) random, (II) on the faults, (III) on the low active points, (IV) on the moderately active points, and (V) on the highly active points in the region. The rate of earthquakes was shown to be related to the total activity field over the region of study. Some universal behaviors of the system are shown to be related to the scheme taken for the initial stimuli. The second part of the paper was devoted to the Multi-fractal analysis, which is exploited to extract the spectrum of the Hurst exponent of time series. The time series for each scheme was analyzed separately by multifractal analysis, for all of which the average Hurst exponent is shown to be lower than 0.5. This is an intrinsic property of anti-correlated time series, for which a large rare event is expected to be followed by a small event. The stimulation of highly active regions (in our (12)    www.nature.com/scientificreports/ study, the points with energies M ≥ 4.0 ), a lowest average Hurst exponent is obtained, meaning that we have the strongest anti-correlated system in this case. An overall phase diagram for the model is sketched for all schemes that were considered in this paper.
A support for our model could be the location of the next earthquake in the study area on the Mosha Fault, which is in agreement with the location of the earthquake occurred on May 7th, 2020 (M W 4.9). There are a lot of degrees of freedom that are ignored in our sandpile-based model, like the stochastic onset of movement Figure 10. Calculation of the statistical function F q using Eq. (13) for the dynamics with avalanches in the whole space. The function of F q vs s display power laws F q (s) ∼ s h(q) , where h(q) depend on q. This feature demonstrates that the time series is a multifractal.  Table: The values of exponents ᾱ , δα , τ S , and τ D for the dynamics with avalanches in the whole space, wherein the δα is defined as the width of f (α) , which is the length of the interval between two successive f (α) = 0.3 . Lower Table: The same for the dynamics with avalanches inside the high-resolution box.  www.nature.com/scientificreports/ (which was considered to be identified by a single threshold in each segment) and the conservation of energy/ stress, which is violated in more realistic models, like the Olami-Feder-Christensen (OFC) model of earthquake 60 . Therefore, in case these simplifications fail, our model is not expected to work properly.